% function poisson_ex1
clear all;
clc
% only in main function and mesh_refine function we declare global
% variables, because the geom modification of geom data cost much when the
% computational mesh is very large.
global V T E TE ET posV posT posE;

% init a few global parameters
d =6;  n_dof = 0;
mesh = 1; caseNum = 3; % (1,2,3,)9,11,12

[cr,desc] = get_auxillary_mat(d);

if mesh == 1
    % prepare mesh 1
    [p,e,t] = initmesh('squareg','Hmax',0.55,'Hgrad',1.99);
    T = t(1:3,:)'; 
    if caseNum==11  % [-2,2]x[-2,2]
        V = [p(1,:)' -p(2,:)']*2;
    else % [0,1]x[0,1]
        V = ([p(1,:)' -p(2,:)']+1)/2;
    end
%     [V,T] = refine_uniform(V,T);
elseif mesh==2
    % prepare mesh 2
    [V,T] = equi_tri_lai(6); % [0,1]x[0,1]
    if caseNum==11  % [-2,2]x[-2,2]
        V = 4*V-2;
    end
end

figure('Position',[100 200 800 450]);subplot(1,2,1);
plot_t(V,T);axis tight; % draw mesh at onece

posV = size(V,1); posT = size(T,1); posE = 0;
areas = tdata;   % init the geom data: E,posE,TE,ET,areas


[dofs,n_dof] = sort_dof_dis(T,posT,d);

% prepare for system matrix
M = mass_bb(dofs,V,T,posT,d,areas);
A = stiff_bb(dofs,V,T,posT,d,areas,desc);
b = M*bnet(dofs,V,T,posT,d,'elliptic_f',caseNum);

% dirichlet boundary condition:
bdr = find(ET(1:posE,2)==0);
[B,bf] = dirichlet_bc_bb(dofs,n_dof,V,E,T,ET,TE,cr,d,bdr,'elliptic_g',caseNum);

 [H,row_idx] = smooth_C1_dis(dofs,V,T,E,posE,ET,TE,d,cr);
%[H,row_idx] = smooth_C1_dis(dofs,V,T,E,posE,ET,TE,d,cr);
x = lagrange22(A,b,[B;H],[bf;zeros(size(H,1),1)],1e-6,15,1e-15);

% eval at given points
[z,tris,pts] = heval(dofs,V,T,posT,x,d,make_mesh(11),11);
u = feval('elliptic_u',pts(:,1),pts(:,2),caseNum);

% visualize the fe solution
subplot(1,2,2);trisurf(tris,pts(:,1),pts(:,2),abs(z-u));